module AssemblyLinearSystem
    implicit none

	integer, parameter :: double=8

    type, abstract :: LinearSystemAssembler

    contains
    	procedure, public :: Adjusted       !< Check if all the prerequisites are on order
    	procedure, public :: DropMatrix     !< Return the extended Matrix of the System
    end type

	abstract interface
		pure function Adjusted(this) result (boolOk)
			import
			class(LinearSystemAssembler), intent(in) :: this
			logical :: boolOk
		end function

		pure function DropMatrix(this, geometry, boundaryData) result (doubleM)
			import
			class(LinearSystemAssembler), intent(in) :: this
			real(kind=double), dimension(:,:), intent(in) :: geometry
			real(kind=double), dimension(:,:), intent(in) :: boundaryData
			real(kind=double), dimension(:,:), allocatable :: doubleM
		end function
	end interface

contains

end module


module NonlinearTransientTemperatureLinearConductivityModel
	use AssemblyLinearSystem
	use iso_varying_string
	implicit none

	integer(kind=1), parameter :: maxPropertiesCount = 4


	!> \brief NonLinear Heat Transfer Temperature with a Linear
	!> Temperature-based Thermal Conductivity Model
	!!
	!! This is a Class for the Assembly of the Element Matrix
	type, public, extends(LinearSystemAssembler) :: NLTLinearConductivity
		private
		type(iso_varying_string) :: MaterialName
		logical, dimension(maxPropertiesCount) :: SettingsSystem
		real(kind=8) :: rho
		real(kind=8) :: cond_LinearCoef
		real(kind=8) :: cond_AngularCoef
		real(kind=8) :: timeStep
	contains
    	procedure, public :: Adjusted       !< Check if all the prerequisites are on order
    	procedure, public :: DropMatrix     !< Return the extended Matrix of the System
		procedure, public :: SetConductivity_LinearCoeficient
		procedure, public :: SetConductivity_AngularCoeficient
		procedure, public :: SetTimeStep
		procedure, public :: SetMassDensity
	end type

contains

	!> \brief Implementation of the Adjusted
	!!
	!! \return boolOk - The verification result
	pure function Adjusted (this) result (boolOk)
		implicit none
		type(NLTLinearConductivity), intent(in) :: this
		logical :: boolOk = .true.
		integer(kind=1) :: i

		do i = 1, maxPropertiesCount
			if (.not. SettingsSystem(i) ) then
				boolOk = .false.
				return
			end if
		end do
	end function


	pure function DropMatrix(this, geometry, boundaryData) result (doubleM)
		implicit none
		type(NLTLinearConductivity), intent(in) :: this
		real(kind=8), dimension(:,:), intent(in) :: geometry
		real(kind=8), dimension(:,:), intent(in) :: boundaryData

		real(kind=8), dimension(:,:), allocatable :: doubleM


	end function





    function GetGradientTemperatureMatrix( Kappa, Rho, Cp, StabilizationFlux, StabilizationTemperature, Temperatures, &
            HeatFluxes_X, HeatFluxes_Y, HeatFluxes_Z) result (Matrix)
        integer, parameter :: Self = 1
        implicit none
        real(kind=8), intent(in) :: kappa
        real(kind=8), intent(in) :: Rho
        real(kind=8), intent(in) :: Cp
        real(kind=8), intent(in) :: StabilizationFlux
        real(kind=8), intent(in) :: StabilizationTemperature
        real(kind=8), dimension(:,:), intent(in) :: Temperatures
        real(kind=8), dimension(:,:), intent(in) :: HeatFluxes_X
        real(kind=8), dimension(:,:), intent(in) :: HeatFluxes_Y
        real(kind=8), dimension(:,:), optional, intent(in) :: HeatFluxes_Z

    end function

    subroutine AssemblyAssociatedHeatFluxTempLinearRectangle(El, diffusion, capacitance, density)
        implicit none

        type(Element), intent(out) :: El
        real(kind=8), intent(in) :: diffusion, capacitance, density

        integer :: i, j

    end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Common Use on the manipulation of TriangularElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    !> \brief Base functions on Triangular Linear coordinate $( \eta_1, \eta_2, \eta_3 )$ system
    !!
    !! \param i - the index of the base function
    !! \param eta - double precision 3d vector of the point on which the base is calculated
    !! \return v - double precision scalar calculated value
    !!
    pure function TriangularPhi(i, eta ) result (v)
        implicit none
        integer, intent(in) :: i
        real(kind=8), dimension(3), intent(in) :: eta
        real(kind=8) :: v
		if(size(eta) /= 3 ) return
        select case(i)
            case(1)
                v = eta(1)
            case(2)
                v = eta(2)
            case(3)
                v = eta(3)
            case default
                v = 0.0d0
        end select
    end function


    !> \brief Base functions on single variable, or linear, system
    !!
    !! \param i - the index of the base function
    !! \param eta - double precision scalar point on which the function is calculated
    !! \return v - double precision scalar calculated value
    !!
    pure function LinearPhi(i, eta) result (v)
        implicit none
        integer, intent(in) :: i
        real(kind=8), intent(in) :: eta
        real(kind=8) :: v

        select case(i)
            case(1)
                v = 0.5d0 * (1.0d0 - eta)
            case(2)
                v = 0.5d0 * (1.0d0 + eta)
            case default
                v = 0.0d0
        end select
    end function

    !> \brief Returns the Jacobian over the 1-D Boundary, or side of the Triangle, of integrals performed
    !!
    !! \param X - the Nodes Coordinates
    !! \param indexes - 3D integer vector comprising the indexex on the rotation of sides on which calculate the Jacobian
    !! \return v - double precision value of the Jacobian
    !!
    pure function JacobianLength(X, indexes) result (v)
		implicit none
        real(kind=8), dimension(nNodes,probDim), intent(in) :: X
        integer, dimension(3), intent(in) :: indexes
        real(kind=8) :: v
        v = 0.5d0*( (X(indexes(3),1) - X(indexes(2),1)**2.0d0 + (X(indexes(3),2) - X(indexes(2),2)**2.0d0)**0.5d0
    end function


    !> \brief Returns a analytically calculated Matrix which comprises of the Temperature contribution within the element area
    !!
    !! \param this - the Element structure
    !! \param X - the nodes coordinates
    !! \param orientation - the Axis representing the Allocation of the Matrix
    !! \return local - The Analytically calculated matrix
    !!
	pure function MatrixH_l( this, X, orientation ) result (local)
		implicit none
		type(TriangleElement), intent(in) :: this
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer, intent(in) :: orientation

		real(kind=8), dimension(nNodes,nNodes) :: local
!// \todo (ronaldo#1#): É Necessario adequar para a condutividade térmica em função da temperatura.

		select case(orientation)
			case(XAxis)
				local (1,1) = X(2,2) - X(3,2) !< $y_{23}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(3,2) - X(1,2) !< $y_{31}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(1,2) - X(2,2) !< $y_{12}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)
			case(YAxis)
				local (1,1) = X(3,1) - X(2,1) !< $x_{32}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(1,1) - X(3,1) !< $x_{13}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(2,1) - X(1,1) !< $x_{21}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)
			case default
				local (1,1) = X(2,2) - X(3,2) !< $y_{23}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(3,2) - X(1,2) !< $y_{31}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(1,2) - X(2,2) !< $y_{12}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)

		end select
		local = local/6.0d0
	end function


    !> \brief Generate a element of a Matrix $E_{lr,i}$, or a Matrix $E_{lr,i,B}$
    !!  Element boundary thermal conductivity contribution on side $i$
    !!  To gerenate a $E_{lr,i,B}$ matrix usa a negative Stabilization parameter
    !! \param this -> the element
    !! \param l, r -> ranges from x to y, each
    !! \param i -> the boundary
    !! \param k, m -> Matrix internal indexes
    !! \param ThermCond -> Thermal Conductivity. Function which defines the Thermal Conductivity
    !! \param Stabilization -> Stabilization Parameter
    !! \return cell - double precision scalar
    !!
	pure function MatrixE_lr_i_km( this, l, r, i, k, m, X ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = 3
		integer, intent(in) :: l, r, i, k, m
        real(kind=8), dimension(:,:), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface

        real(kind=8) :: cell
		integer :: g !< integration points index
		integer :: a !< Thermal Conductivity interpolation Index

		real(kind=8), dimension(IntPoints) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		integer, dimension(3) :: vertex

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (a=2:3)             !< Points not regarding the vertex pivoting the situation
			forall (g=1:IntPoints) !< Warp through Integration Points
				cell = cell + ThermCond(this%DOF(Temperature_Initial+a-1))*&
							LinearPhi(vertex(a), IntegrationPoints(g)) *&
							LinearPhi(k, IntegrationPoints(g)) *&
							LinearPhi(m, IntegrationPoints(g)) *&
							IntegrationWeights(g)
			end forall
		end forall
		cell = this%OuterVector(i,l) * this%OuterVector(i,r) * JacobianLength(X, vertex) * Stabilization * cell
	end function


    !> \brief Returns a Component cell of the Matrix $H_{l,i}$ Calculated over the side $i$ of the element
    !!
    !! \param this - the Element Structure
    !! \param l - the orientation Axis of the Calculated Matrix
    !! \param i - the side on which this calculations are being performed
    !! \param k - the local matrix $H_{l,i}$ line of the component
    !! \param m - the local matrix $H_{l,i}$ column of the component
    !! \param X - the Nodes components
    !! \param ThermCond - the function defining the Thermal Conductivity at node $n$ depending, only, on temperature $T_n$
    !! \return cell - double precision scalar
    !!
	pure function MatrixH_l_i_km( this, l, i, k, m, X, ThermCond ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: IntPoints = 1
        type(TriangleElement), intent(in) :: this
		integer, intent(in) :: l, i, k, m
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface

        real(kind=8) :: cell
		integer :: g !< integration points index
		integer :: a !< Thermal Conductivity interpolation Index

		real(kind=8), dimension(IntPoints,3) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (a=2:3) !< Points not regarding the vertex pivoting the situation
			forall (g=1:IntPoints) !< Warp through Integration Points
				cell = cell + ThermCond(this%DOF(Temperature_Initial+a-1))*&
							LinearPhi(vertex(a), IntegrationPoints(g)) *&
							LinearPhi(k, IntegrationPoints(g)) *&
							LinearPhi(m, IntegrationPoints(g)) *&
							IntegrationWeights(g)
			end forall
		end forall
		cell = cell  * JacobianLength(X, vertex) * this%OuterVector(i,l)
	end function


    !> \brief Returns a component of the Matrix $J_{l,i}$ calculated over the side $i$ of the element
    !!
    !! \param this - the Element Structure
    !! \param l - the orientation Axis of the Calculated Matrix
    !! \param i - the side on which this calculations are being performed
    !! \param k - the local matrix $J_{l,i}$ line of the component
    !! \param m - the local matrix $J_{l,i}$ column of the component
    !! \param X - the Nodes components
    !! \return cell - double precision scalar
    !!
	function MatrixJ_l_i_km( this, l, i, k, m, X ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = 3
		integer, parameter :: nDimension = 2
        integer, parameter :: IntPoints = 1
        type(TriangleElement), intent(inout) :: this
		integer, intent(in) :: l, i, k, m
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element

        real(kind=8) :: cell
		integer :: g !< integration points index

		real(kind=8), dimension(IntPoints) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		integer, dimension(3) :: vertex

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (g=1:IntPoints) !< Warp through Integration Points
			cell = cell + LinearPhi(k, IntegrationPoints(g)) *&
						LinearPhi(m, IntegrationPoints(g)) *&
						IntegrationWeights(g)
		end forall
		cell = cell * JacobianLength(X, vertex) * this%OuterVector(i,l)
	end function


    !> \brief Analitically calculated Matrix component of the block-matrices $G_{T,i}$ and $G_{T,i,B}$
    !!
    !! \param i - the side on which the calculation is being performed
    !! \param X - The Nodes Coordinates
    !! \return C - The Matrix itself
    !!
   pure function MatrixC_i( i, X ) result (C)
        implicit none
        integer, intent(in) :: i
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

        real(kind=8), dimension(nNodes,nNodes) :: C

        integer, dimension(nNodes) :: vertex
        integer :: a, b

        C = 0.0d0
        vertex = rotation(i)
		C(vertex(2),vertex(2)) = 2.0d0
		C(vertex(3),vertex(3)) = 2.0d0
		C(vertex(2),vertex(3)) = 1.0d0
		C(vertex(3),vertex(2)) = 1.0d0
        C = C * ( &
					(X(vertex(2),1)-X(vertex(3),1))**2.0d0 + &
					(X(vertex(2),2)-X(vertex(3),2))**2.0d0 &
					)**0.5
    end function




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Use on the calculation of the Stiffness Matrix !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !> \brief Auxiliary function S, which gives the ponderation value for the calculations of the fixed-form, or analytically calculated matrices on the triangle
    !!
    !! \param X - the Nodes coordinates
    !! \return local - the value
    !!
	pure function S (X) result (local)
		implicit none
		real(kind=8), dimension(nNodes, probDim), intent(in) :: X
		real :: local
		local = ( (X(1,1)-X(3,1))*(X(2,2)-X(3,2)) - (X(2,1)-X(3,1))*(X(1,2)-X(3,2)) )/2.0d0
	end function


    !> \brief Returns the Elementary matrix E, refered to the block-diagonal of the Stiffness local Matrix
    !!
    !! \param X - the Nodes coordinates
    !! \return local - The Entire Matrix
    !!
	pure function Matrix_E(X) result (local)
		implicit none
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		real(kind=8), dimension(nNodes,nNodes) :: local
		local = 1.0d0
		local(1,1) = 2.0d0; local(2,2) = 2.0d0; local(3,3) = 2.0d0
		local = local * S(X)/12.0d0
	end function


    !> \brief Set the E central matrices of the local matrix
    !! Matrix of the inner element contribution relative to the heat flux
    !! The E matrices are dependent only of a determinant obtained from the Vertices Coordinates
    !! The E matrices compond the diagonal blocks of the first assembly
    !! \param this - The Triangle Element itself
    !! \param S - The determinant parameter
    !!
    subroutine SetMatricesE( Matrix, X )
        implicit none
		real(kind=8), dimension(MatrixSize,MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l, i
		integer, dimension(nNodes) :: idx

		forall(l = 1:probDim, i=i:nNodes)
			forall(i=1:nNodes) !< Set Range for Matrices E
				idx(i) = nNodes*(l-1)+i+1
			end forall
			idx = SetRange(l)
			Matrix(idx,idx) = Matrix(idx,idx) + Matrix_E(X)
		end forall
    end subroutine


    !> \brief Set the $H_k$ matrix of the local matrix, $k = x,y$
    !! Matrix of the inner element contribution relative to Temperature
    !! \param this - The Triangular Element itself
    !! \param X - The Coordinates of the nodes
    !!
	subroutine SetMatricesH_l( this, Matrix, X )
		implicit none
		type(TriangleElement), intent(inout) :: this
		real(kind=8), dimension(MatrixSize, MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer :: l

		forall( l = 1:probDim)
			Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) = Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) + MatrixH_l( this, X, l)
		end forall

	end subroutine


    !> \brief Set the $J_k$ matrix of the local matrix, $k = x,y$
    !! Matrix of the inner element contribution relative to Temperature
    !! \param this - The Triangular Element itself
    !! \param X - The Coordinates of the nodes
    !!
	subroutine SetMatricesJ_l( this, Matrix, X )
		implicit none
		type(TriangleElement), intent(inout) :: this
		real(kind=8), dimension(MatrixSize, MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer :: l

		forall( l = 1:probDim)
			Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) = Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) + MatrixH_l( this, X, l)
		end forall

	end subroutine


	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Face components of the Stiffness Matrix!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !> \brief Set the $E_{lr,i}$ matrices comprising the element contribuitions to the Stiffness matrix on every side of the element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param ThermCond - The Thermal Conductivity function which determines the isotropic thermal conductivity coeficient
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
    subroutine SetMatricesE_lr_i( this, Matrix, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize,MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		integer :: l, r  !< l, r = x, y, z (1, 2, 3)
        integer :: i !< side index
        integer :: k, m !< local base index
		integer :: l_m, r_m  !< local matrix real position

        forall(i = 1:nSides)  !< Warp through  Triangle1's Sides
			forall (l=1:nDimension) !< Warp through Local dimensions
				forall (r=1:nDimension) !< Warp through Local dimensions
					forall (k=1:nInterpolationPoints) !< Warp through  local basis
						forall (m=1:nInterpolationPoints) !< Warp throug local basis
							l_m = (l-1)*nNodes+k  !< Matrix Index of Element Matrix
							r_m = (r-1)*nNodes+m  !< Matrix Index of Element Matrix
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + &
									MatrixE_lr_i_km(this, l, r, i, k, m, X, ThermCond, Stabilization) *&
									Stabilization
						end forall
					end forall
				end forall
			end forall
        end forall
    end subroutine


    !> \brief Sets the $H_{l,i}$ matrices comprising the Contributions of Temperature on each side of the Element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param ThermCond - The Thermal Conductivity function which determines the isotropic thermal conductivity coeficient
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
	subroutine SetMatricesH_l_i( this, Matrix, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize, MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + ( -(0.5d0 + Stabilization)) * &
									MatrixH_l_i_km(this, l, i, k, m, X, ThermCond)
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the $J_{l,i}$ Temperatures comprising the Heat Fluxed contribution to the temperature on each side of the element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
	subroutine SetMatricesJ_l_i( this, Matrix, X, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize, MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + ( -(0.5d0 - Stabilization)) * &
									MatrixJ_l_i_km(this, l, i, k, m, X, Stabilization)
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the $G_{T,i}$ matrices of the Stiffness Matrix comprising the Temperature Contributions to the Temperature itself on each side
    !!
    !! \param Matrix - The matrix which the calculations will be recorded
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
    subroutine SetMatrixG_T_i( Matrix, X, Stabilization )
        implicit none
        integer, parameter :: nSides = nNodes
		real(kind=8), dimension(MatrixSize,MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l

		forall(l = 1:nSides)
			Matrix(idx,idx) = Matrix + Stabilization * MatrixC_i(l, X)
		end forall
    end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!! Use on the calculations of the Neighbour dominated vector of the Element !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



    !> \brief Set the $E_lr_i_B$ comprising the Neighbours contribution to the independent constant vector due to the temperature effect on the flux
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !! \param ThermCond - The Thermal Conductivity function for update the values
    !!
    subroutine SetMatricesE_lr_B_i( this, Vector, Neighbours, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		integer :: l, r  !< l, r = x, y, z (1, 2, 3)
        integer :: i !< side index
        integer :: k, m !< local base index
		integer :: l_m, r_m  !< local matrix real position

        forall(i = 1:nSides)  !< Warp through  Triangle1's Sides
			forall (l=1:nDimension) !< Warp through Local dimensions
				forall (r=1:nDimension) !< Warp through Local dimensions
					forall (k=1:nInterpolationPoints) !< Warp through  local basis
						forall (m=1:nInterpolationPoints) !< Warp throug local basis
							l_m = (l-1)*nNodes+k  !< Matrix Index of Element Matrix
							r_m = (r-1)*nNodes+m  !< Matrix Index of Element Matrix
							Vector(l_m) = Vector(l_m) - &
									MatrixE_lr_i_km(this, l, r, i, k, m, X, ThermCond, Stabilization) *&
									Stabilization * Neighbours(i)%DoF(r_m, Dof_Commit) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
        end forall
    end subroutine


    !> \brief Sets the matrices $H_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the temperature effect on the flux
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Matrix - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
	subroutine SetMatricesH_l_B_i( this, Vector, Neighbours, X, Stabilization ) !< Linear Temperature Conductivity
        implicit none
        type(NLTLinearCondutivity), intent(inout) :: this          !< Nonlinear Transient Linear Temperature-dependent Conductivity
        real(kind=8), dimension(:,:), intent(inout) :: Matrix      !< Destination Matrix
        real(kind=8), dimension(:,:), intent(in) :: Neighbours     !< Values of DoF in the neighbours
        real(kind=8), dimension(:,:), intent(in) :: X              !< Values of Nodes Positions
        real(kind=8), intent(in) :: Stabilization                  !< Stabilization Parameter

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m
		integer :: nSides, nDimension, nInterpolationPoints

		nSides = dim(X,1)                         !< Number of Sides for
		nDimension = dim(X,2)                     !< Dimension of the problem
		nInterpolationPoints = dim(Neighbours, 1) !< Interpolation 'Nodal' Points

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Vector(l_m) = Matrix(l_m, r_m) + ( -(0.5d0 - Stabilization)) * &
									MatrixH_l_i_km(this, l, i, k, m, X, ThermCond) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the matrices $J_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the flux effect on the temperature
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
	subroutine SetMatricesJ_l_B_i( this, Vector, Neighbours, X, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Vector(l_m) = Vector(l_m) + ( -(0.5d0 + Stabilization)) * &
									MatrixJ_l_i_km(this, l, i, k, m, X) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the matrices $G_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the flux effect on the temperature
    !!
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
    subroutine SetMatrixG_T_B_i( Vector, Neighbours, X, Stabilization )
        implicit none
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l

		forall(l = 1:nSides)
			Matrix(idx,idx) = Matrix - Stabilization * MatrixC_i(l, X) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
		end forall
    end subroutine





end module
